districtcovariates="pct_hisp_2005:factor(year) + pct_black_2005:factor(year) + pct_male_2005:factor(year) + pct_white_2005:factor(year) + enroll_2005:factor(year) + pct_ell_2005:factor(year) + pct_speced_2005:factor(year) + pct_FRPM_2005:factor(year) + missing_pct_ell_2005 + missing_pct_speced_2005 + missing_pct_FRPM_2005 "



acs="pctemploy2005:factor(year)  + pctsinglemother2005:factor(year) + pctmanufac2005:factor(year) + pctutility2005:factor(year) + pctbachelorhigher2005:factor(year)"
# import acs
#############################
# Table A.4 Panel B: 
distance=60
filetoimport=paste("Derived Data/dfbartik2005_",distance,"kmanalysis.csv",sep="")

dfbartik=fread(filetoimport)
dfbartik=merge(dfbartik,tempdf,by=c("leaid","year"),all.x=TRUE)

dfbartik=dfbartik %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort,pctutility2005,pctemploy2005, pctsinglemother2005,pctmanufac2005, pctutility2005, pctbachelorhigher2005, pct_hisp_2005,pct_male_2005,enroll_2005,pct_ell_2005,pct_speced_2005,pct_FRPM_2005, pct_black_2005, pct_white_2005)

sample_taba4_b=dfbartik[dfbartik$testmonth=="May",]
# instrument: 
instname=paste("bartikleave_",distance,sep="")
# weights:
weightval=paste("totalprod_fuel",distance,sep="")


#  district covariates
bartik_dist=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))

resultbartik_dist=felm(bartik_dist,sample_taba4_b,weights=sample_taba4_b[[weightval]],na.action=na.omit)


bartik_sort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ",districtcovariates, "+ ", acs ," |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~",instname,")|leaid",sep="" ))
resultbartik_sort=felm(bartik_sort,sample_taba4_b,weights=sample_taba4_b[[weightval]],na.action=na.omit)


# district + acs +weather long 
bartik_weathavg=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates, " + ", acs ," + ", tempvarsavg ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
resultbartik_weathavg=felm(bartik_weathavg,sample_taba4_b,weights=sample_taba4_b[[weightval]],na.action=na.omit)


# export 
fileout= paste("output/TableA4_PanelB.tex",sep="")

stargazer(resultbartik_dist,resultbartik_sort,resultbartik_weathavg, keep="avgpm25_9",digits=4, add.lines=list(c("First-Stage F Stat",resultbartik_dist$stage1$iv1fstat$avgpm25_9[5], resultbartik_sort$stage1$iv1fstat$avgpm25_9[5], resultbartik_weathavg$stage1$iv1fstat$avgpm25_9[5]),
                                                                                                               c("Student Covariates", "Yes", "Yes", "Yes"),
                                                                                                               c("Local Economic Controls","No", "Yes", "Yes"),
                                                                                                               c("Weather Controls",   "No", "No", "Yes")
),
out=fileout )
###############################################
# Table A.4 Panel C

# Table A.4 Panel C: 
distance=40
filetoimport=paste("Derived Data/dfbartik2005_",distance,"kmanalysis_annual.csv",sep="")

dfbartik=fread(filetoimport)
dfbartik=merge(dfbartik,tempdf,by=c("leaid","year"),all.x=TRUE)

dfbartik=dfbartik %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort,pctutility2005,pctemploy2005, pctsinglemother2005,pctmanufac2005, pctutility2005, pctbachelorhigher2005, pct_hisp_2005,pct_male_2005,enroll_2005,pct_ell_2005,pct_speced_2005,pct_FRPM_2005, pct_black_2005, pct_white_2005)

sample_taba4_c=dfbartik
# instrument: 
instname=paste("bartikleave_",distance,sep="")
# weights:
weightval=paste("totalprod_fuel",distance,sep="")


#  district covariates
bartik_dist=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))

resultbartik_dist=felm(bartik_dist,sample_taba4_c,weights=sample_taba4_c[[weightval]],na.action=na.omit)


bartik_sort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ",districtcovariates, "+ ", acs ," |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~",instname,")|leaid",sep="" ))
resultbartik_sort=felm(bartik_sort,sample_taba4_c,weights=sample_taba4_c[[weightval]],na.action=na.omit)


# district + acs +weather long 
bartik_weathavg=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates, " + ", acs ," + ", tempvarsavg ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
resultbartik_weathavg=felm(bartik_weathavg,sample_taba4_c,weights=sample_taba4_c[[weightval]],na.action=na.omit)


# export 
fileout= paste("output/TableA4_PanelC.tex",sep="")

stargazer(resultbartik_dist,resultbartik_sort,resultbartik_weathavg, keep="avgpm25_9",digits=4, add.lines=list(c("First-Stage F Stat",resultbartik_dist$stage1$iv1fstat$avgpm25_9[5], resultbartik_sort$stage1$iv1fstat$avgpm25_9[5], resultbartik_weathavg$stage1$iv1fstat$avgpm25_9[5]),
                                                                                                               c("Student Covariates", "Yes", "Yes", "Yes"),
                                                                                                               c("Local Economic Controls","No", "Yes", "Yes"),
                                                                                                               c("Weather Controls",   "No", "No", "Yes")
),
out=fileout )


